† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant Nos. 11622437, 61674171, 11804247, and 11974422), the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB30000000), Key Research Program of Frontier Sciences, Chinese Academy of Sciences (B.L, W.Z.), the Fundamental Research Funds for the Central Universities, China, and the Research Funds of Renmin University of China [Grant Nos. 16XNLQ01 and No. 19XNQ025 (W.J.)].
Point and line defects are of vital importance to the physical and chemical properties of certain two-dimensional (2D) materials. Although electron beams have been demonstrated to be capable of creating single- and multi-atom defects in 2D materials, the products are often random and difficult to predict without theoretical inputs. In this study, the thermal motion of atoms and electron incident angle were additionally considered to study the vacancy evolution in a black phosphorus (BP) monolayer by using an improved first-principles molecular dynamics method. The P atoms in monolayer BP tend to be struck away one by one under an electron beam within the displacement threshold energy range of 8.55–8.79 eV, which ultimately induces the formation of a zigzag-like chain vacancy. The chain vacancy is a thermodynamically metastable state and is difficult to obtain by conventional synthesis methods because the vacancy formation energy of 0.79 eV/edge atom is higher than the typical energy in monolayer BP. Covalent-like quasi-bonds and a charge density wave are formed along the chain vacancy, exhibiting rich electronic properties. This work proposes a theoretical protocol for simulating a complete elastic collision process of electron beams with 2D layers and will facilitate the establishment of detailed theoretical guidelines for experiments on 2D material etching using focused high-energy electron beams.
Since graphene monolayers were first exfoliated,[1–5] two-dimensional (2D) materials have received extensive attention[6–11] owing to their unique physical and chemical properties as well as potential applications.[12–19] Defects are inevitable in 2D materials and may play key roles in tailoring their physical and chemical properties.[20–26] For example, vacancies and domain boundaries, as two typical types of defects in transition metal dichalcogenide (TMDC) monolayers, often show specific local electronic[21] and magnetic structures[21,23] and even localized excitons.[25] However, it has long been a challenge to introduce and modify functional defects in a feasible and controllable manner.[27–30]
Focused electron beams have emerged as a possible means of manipulating individual atoms on the atomic scale, benefitting from the development of scanning transition electron microscopy (STEM). Various defects, even those stabilized only under extremely non-equilibrium conditions, such as multiple vacancies in graphene[31,32] and anti-site defects in MoS2,[21,33] have been successfully obtained using focused high-energy electron beams (FHEEBs) with kinetic energies greater than 60 keV.[27] The stable non-equilibrium state is an inherent benefit of material preparation using FHEEBs rather than conventional synthesis methods.[34–39] FHEEBs are also capable of creating much larger defects from intact layers, e.g., boundaries in TMDCs[39] and narrow nanoribbons of graphene.[36] However, the vacancies created in 2D materials, e.g., graphene and TMDCs,[22,31] are usually randomly distributed, and only the statistics method could be employed to perform analysis.[21,34,36] In other words, their formation process is complicated, and it is very expensive to capture the dynamic trajectories experimentally. In light of these challenges, insights from theoretical simulations have been of paramount importance to the field of atomic manipulation using FHEEBs.
A 2D material of great interest, black phosphorus (BP) has recently been found to exhibit excellent mechanical and electronic properties,[40–45] which could also be manipulated by FHEEBs.[46] A previous work revealed the dynamic stabilities of single-atom defects and edge atoms under FHEEBs.[47] However, the development of these defects and edges was not considered in detail, although this development determines the final type of the defect formed. In this paper, we focus on the evolution of a defect from a single-atom vacancy under a FHEEB in monolayer BP and predicted the final formation of a chain vacancy, composed of two parallel zigzag edges along the (100) direction. During the evolution, P atoms in the same zigzag chain along the (100) direction tend to be struck away one by one because the displacement threshold energies (Td) of these atoms are 8.55–8.79 eV, always lower than those along other directions during the vacancy evolution. The zigzag-like chain vacancy exhibits rich electronic properties. Firstly, there are clear antibonding- and bonding-like states along the chain vacancy, indicating that the interaction between the two edges of the vacancy is a new kind of covalent-like quasi-bonding (CLQB). Secondly, structural ups and downs are formed with an obvious bandgap of 0.13 eV opened, which is a typical characteristic of a charge density wave (CDW). This work provides a more comprehensive protocol for modelling the elastic collision process of high-energy electrons with 2D layers and is expected to inspire future works that will provide detailed theoretical guidance for experiments in the field of 2D material etching by FHEEBs.
Our density functional theory (DFT) calculations were performed using a generalized gradient approximation for the exchange–correlation potential, the projector augmented wave method,[48,49] and a plane-wave basis set implemented in the Vienna ab-initio simulation package (VASP). Van der Waals interactions were considered at the vdW-DF level with the optB86b functional for exchange (optB86b-vdW),[50] which was found to be suitable for modelling the structural properties of BP.[40] The kinetic energy cut-off for the plane-wave basis set was set to 500 eV for the structural relaxation calculations. A 2 × 2 × 1 k-mesh was adopted to sample the first Brillion zone of a 7 × 5 supercell of monolayer BP. The shape and volume of the primitive cell were fully optimized, and all of the atoms in the supercell were allowed to relax until the residual force per atom was less than 0.01 eV/Å. The formation energy of the BP edges was derived from EForm = (EEdge − NPμP)/NP, where EEdge is the total energy of the edge and μP represents the chemical potential of one P atom in monolayer BP.
In previous works, an ab-initio molecular dynamics (AIMD) method was proposed to simulate electron beam irradiation.[32,35,47,51] The effects of elastic collisions were assessed by evaluating the displacement threshold energy Td, which was defined as the minimum initial kinetic energy needed to displace an atom (referred to as a target atom in this paper) from its original site without it recombining with the vacancy. Although a finite temperature environment was previously considered by introducing the Debye model, the random thermal motions were not actually simulated.[32]
In this study, we considered the random thermal motions by introducing a canonical ensemble process at 300 K before performing simulations. To obtain Td, AIMD was used to simulate the structural response to incoming momentum transferred from the incident high-energy electron beams. The kinetic energy cut-off was set to 400 eV in the AIMD simulations.
The simulation process consisted of two steps. Firstly, a canonical ensemble molecular dynamics simulation was performed at 300 K for 4 ps with a time step of 1 fs to simulate a thermal vibration environment. To reduce the effects of velocity fluctuations, two extreme steps were conducted in the last 2 ps to continue the simulation, where the target P atom had the highest velocities along and against the normal vector of the BP monolayers (which could be defined as the z-direction). Secondly, in both cases, we introduced an additional initial velocity along the −z direction toward the target P atom to simulate the instantaneous momentum transfer from the downward electron beam. We then performed a microcanonical ensemble simulation for 0.8–2 ps with a time step of 0.5 fs and traced the trajectory of the target P atom to judge whether the atom was out of the lattice structure. Two criteria were used here: the displacement of the target atom and the displacement times the velocity of the target atom. The first criterion was set to 5 Å, and the second was estimated as 2 Å × 0.05 Å/fs because the interaction was too weak to attract the target P atom back to the original site in these cases. A binary search method was used to obtain the range of the transferred kinetic energy until the velocity resolution was less than 0.001 Å/fs (Td ∼ 0.1 eV). Finally, Td was calculated from the average of the threshold velocities in the two extreme cases.
A case with oblique electron beam incidence was also considered by introducing additional velocities tilted at 10° to the target atoms in the pristine BP monolayer. Oblique incidence did not effectively reduce the displacement threshold energies of the P1T1 atoms, which also had difficulty escaping. For P1D1, Td increased by 0.3–1.0 eV (tilted along four different directions (±1, 0, 0) and (0, ±1, 0)) compared to that in the vertical incidence case. To unify the standards, we evaluated the minimum Td in the above investigations, so only the vertical cases were considered.
The entire simulation process involved extensive labor. A control program, called the automatic beam effect simulation tool (aBEST), was thus written to implement the above process automatically by calling VASP within Python language. The aBEST and example files are accessible through https://gitee.com/jigroupruc/aBEST.
Performances of the optB86b-vdW, PBE-D3, and PBE functionals were examined by calculating Td for P1D2 and P2D1 of configuration V1P. The inclusion of van der Waals corrections reduces the estimated Td, i.e., from 6.68 eV and 8.55 eV of PBE to 6.57 eV and 7.10 eV of PBE-D3 and 5.88 eV and 7.42 eV of optB86b-vdW for atoms P1D2 and P2D1. However, their relative differences between atoms P1D2 and P2D1 are rather comparable, especially for those of PBE and optB86b-vdW. In light of this, the use of PBE, instead of optB86b-vdW, in DFT-MD simulations does not qualitatively change the conclusion of relative Td energies of different atoms, but saves the computational cost by roughly three times.
The electron beam threshold kinetic energy here is defined as the minimum electron beam kinetic energy required to knock out the target atom. In other words, the cross-section of the electron beam for the target atom becomes positive just at the threshold kinetic energy. The cross-section at a selected Td can be derived using the McKinley–Feshbach formalism[52] as
There are two kinds of P atoms in pristine monolayer BP when electron beams are incident from the top to bottom because of the symmetry of the geometry, as shown using different colors in Fig.
The Td was calculated for two typical atoms (P1T1, P1D1) to compare their dynamic stabilities under the FHEEB. It is very difficult (Td > 19 eV) for the P1T1 atom (see Fig.
The chemical environments of the atoms surrounding the single-atom vacancy are very different from those in pristine monolayer BP. Td was calculated for five selected typical atoms to uncover the dynamic stability differences of the surrounding atoms. These five atoms represent three typical directions in which the vacancy could develop, as indicated in Fig.
Continuing the vacancy development, Td and
The reason for this finding is associated with the local structures around the atoms, according to our observation of the atom trajectories. The atoms in the upper sublayers always have difficulty escaping because they would bond to or crash into the atoms in the lower sublayers although the electron beam breaks their bonds with other upper-layer atoms. Meanwhile, for the atoms in the lower sublayers, there are fewer atoms nearby to prevent them from escaping when they are near vacancies than in the pristine lattice. In other words, the closer an atom is to a vacancy, the less likely it is to be bonded with and transfer momentum to the surrounding atoms.
Based on these results, we also predicted that the vacancy would expand along a 1D zigzag chain and eventually become a zigzag chain vacancy (as shown in Fig.
The electronic properties of the zigzag chain vacancy were calculated for in-depth study. The differential charge density (DCD) indicated that the interaction between the two parallel zigzag edges is stronger than the van der Waals interaction because the electron charge significantly decreases near the interfacial P atoms and accumulates between them, as shown in Fig.
A CDW feature was also found in this zigzag chain vacancy. The cross-point of the Fermi level and MB2 was near the midpoint between points Γ and X. A double-periodic lattice was thus tested to look for the CDW according to the Fermi surface nesting theory in a one-dimensional system. The zigzag chain vacancy indeed spontaneously formed structural ups and downs with a tiny energy drop (only 4.3 meV/edge atom), as shown in Figs.
A zigzag edge is obviously not the only kind of edge in monolayer BP. Thus, the geometric structures of several possible edges were relaxed, and their formation energies were calculated to study the thermodynamic stability of our predicted vacancy. Along three crystal directions with low indices, as shown in Figs.
A special zigzag chain vacancy in monolayer BP was predicted by using an FHEEB and knocking away P atoms one by one along a zigzag chain in the lower sublayers. The calculated electronic properties of the chain vacancy showed that there was quasi-bonding between the two edges of the vacancy, and a CDW was also formed along the vacancy. Our findings help improve understanding of quasi-bonding in which covalent-like states can also be half-occupied. The chain vacancy was a dynamically stable but thermodynamically metastable state according to our comparison of the stabilities of five typical edges in monolayer BP. It was inspiring that the electron beam could create a dynamically mostly stable but thermodynamically metastable vacancy, which is difficult to obtain using conventional chemical synthesis methods but easier to achieve using an electron beam, as stated above. This characteristic proves that an FHEEB can create a special environment for defect development. In addition, although the protocol provided in this work offers a route to more comprehensively capture the elastic collision process of defect formation in monolayer BP, the method is also subject to improvement with considering inelastic collisions (or electronic excitation process). Electronic excitation under electron beams is a complicated process which may involve phonon scattering,[56,57] intra- or inter-band transition,[58,59] plasmon,[60] among the others. A few previous methods do be available to simulate those excitation processes, like constrained density functional theory,[61] impulsive 2-state model,[62,63] non-adiabatic molecular dynamics,[64] and time-dependent density functional theory.[65] This work is expected to inspire further works that will implement those exciton modeling methods into the simulation protocol and thus provide detailed theoretical guidance for future experiments in the field of 2D material etching by FHEEBs.
[1] | |
[2] | |
[3] | |
[4] | |
[5] | |
[6] | |
[7] | |
[8] | |
[9] | |
[10] | |
[11] | |
[12] | |
[13] | |
[14] | |
[15] | |
[16] | |
[17] | |
[18] | |
[19] | |
[20] | |
[21] | |
[22] | |
[23] | |
[24] | |
[25] | |
[26] | |
[27] | |
[28] | |
[29] | |
[30] | |
[31] | |
[32] | |
[33] | |
[34] | |
[35] | |
[36] | |
[37] | |
[38] | |
[39] | |
[40] | |
[41] | |
[42] | |
[43] | |
[44] | |
[45] | |
[46] | |
[47] | |
[48] | |
[49] | |
[50] | |
[51] | |
[52] | |
[53] | |
[54] | |
[55] | |
[56] | |
[57] | |
[58] | |
[59] | |
[60] | |
[61] | |
[62] | |
[63] | |
[64] | |
[65] |